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This paper gives an overview of new developments of the least-squares finite element 
method(LSFEM) in fluid dynamics. Special emphasis is placed on the universality of 
LSFEM; the symmetry and positiveness of the algebraic systems obtained from LSFEM; 
the accommodation of LSFEM to equal-order interpolations for incompressible viscous 
flows; and the natural numerical dissipation of LSFEM for convective transport problems 
and high-speed compressible flows. The performance of LSFEM is illustrated by numerical 
examples. 


1. Introduction 

It has been more than two decades since people started applying finite element meth- 
ods to fluid dynamics problemsfl]. Subsequently, a large number of papers and many books 
on finite elements in fluids have been published, and numerous fluid dynamics problems 
have been succesfully solved by finite element methods. Currently, however, only two com- 
mercial general-purpose finite element packages for incompressible flows, i.e. FIDAP and 
FLOTRAN[2], are available. This situation differs from that in solid mechanics. Twenty 
years after the first finite element papers were published, several dozen of finite element 
commercial packages for solid mechanics problems were on the market[3]. Perhaps, one 
reason for the lack of fluid dynamics finite element codes was that there was no unified 
method which could cover a wide range of fluid problems. For example, the classic Galerkin 
method is used for potential flows, the mixed Galerkin method and the penalty method 
are dominant for incompressible viscous flows[4,5], the Taylor-Galerkin method and the 
Petrov-Galerkin method are developed for convective transport problems and compressible 
flow problems [6-8]. Because the principles and structures of these methods are different, 
it is extremely difficult to implement these methods in a general-purpose code. 
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The purpose of this paper is to advocate an unified method - LSFEM for fluid dy- 
namics problems, to provide a conceptual framework for the method, and to review recent 
developments. In addition, comparisons are made with other methods and some represen- 
tative numerical results are presented. 

An outline of this paper follows: In Section 2, a brief review of the literature is given. 
The construction of LSFEM for first-order partial differential equations is presented in 
Section 3. Application of LSFEM for incompressible vicous flows is discussed in Section 

4. Attention is fixed on the accommodation of equal-order interpolations. In Section 

5, LSFEM for a convection problem is described. It is shown that the method natually 
contains an upwinding mechanism to stabilize the solution. In section 6, LSFEM for the 
compressible Euler equations is presented. The capability of LSFEM to capture shocks is 
shown by numerical results. Conclusions are drawn in Section 7. 


2. Review of Literature 

The least-squares method of interest here is based on minimizing the residual in differ- 
ential equations and boundary conditions in a least- squares sense. The use of least-squares 
method for the approximate solution of partial differential equations has a long history [9]. 
The modern theory of least-squares method for the numerical solution of elliptic boundary- 
value problems was started with the papers of Bramble and Schatz[lO,ll], and discussed by 
Varga[l2]. This work uses a finite-dimensional space of approximating functions, similar 
to the spaces used in finite element methods. 

When compared with the classic Galerkin finite element method, LSFEM generally 
leads to more stringent continuity requirements for the trial functions. Lynn and Arya[13] 
and Zienkiewicz et al [14] demonstrated that it is possible to reduce the order of continuity 
requirements at the expense of introducing more unknowns by first transforming the orig- 
inal differential equations into an equivalent system of first-order differential equations. 
Hence, the smoothness requirement for the space of approximate functions would be re- 
duced, and the C° elements would be applicable. Subsequently, LSFEM were applied to 
certain simple flow problems (see Lynn[15], Lynn and Alani [16], and Polk and Lynn[17]). 
Fix and Guzburger [18] also studied LSFEM for systems of first-order equations and results 
were reported for the Tricomi equation. 

A theoretical analysis of LSFEM for the elliptic system of the Petrovsky type was 
developed by Wendland[19], and optimal error estimates were obtained. Petrovsky systems 
are an important subclass of the class of elliptic systems, in which the different equations 
and unknowns appearing in the system have the same “differentiability order” . 

A LSFEM for the Euler equations governing inviscid compressible flows was proposed 
by Fletcher[20]. The feature of his formulation is the designation of groups of variables 
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rather than single variables. Jiang and Chai[21] applied LSFEM to a first-order quasi-linear 
system for compressible potential flow. 

Nguyen and Reynen[22] presented a space-time LSFEM for the advection-diffusion 
problems. Their numerical results show that the use of upwinding techniques or the Taylor- 
Galerkin approach in finite elements turns out to be unnecessary when the least-squares 
weak formulation is extended into the time domain using standard shape functions. 

An error estimate for the least-squares finite element solution of Cauchy-Riemann 
type equations was given by Fix and Rose[23]. 

A more general mathematical theory of LSFEM for the approximate solution of an 
elliptic system of Agmon-Douglis-Nirenberg(ADN) type was presented by Aziz, Kellogg 
and Stephens[24]. The method minimizes a least squares functional that consists of a 
weighted sum of the residuals occuring in the equations and the boundary conditions of 
the system. The weights are determined by the indices that enter into the definition of the 
ADN boundary-value problem. 

It is shown that for the Cauchy-Riemann equations the optimal weighting factor is 
unity [25, 26], The LSFEM for incompressible and compressible flow analysis is implemented 
in conjunction with an element-by-element preconditioned conjugate gradient method and 
an adaptive refinement strategy[25, 27-29]. A new preconditioner based on the H 1 gradient 
was formulated[25,26]. 

The recent advance of application of LSFEM in fluid dynamics will be discussed in 
the following sections. 


3. Construction of LSFEM for First-Order Systems 

Almost all problems arising in fluid dynamics, solid mechanics, heat transfer, electro- 
magnetics and other mathematical physics can be recast in the form of first-order systems. 
For convenience, here we shall restrict our treatment initially to steady-state linear prob- 
lems. It is straigtforward to extend the treatment to quasi-linear problems. 

At first we consider the boundary-value problem: 


Lu — f in fl 

(3.1o) 

Btt = g on T 

(3.16) 

where L is a first-order partial differential operator: 


T V 4 * du * 

L u = } A b A u 

T-* OX; 

(3.2) 
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n E tft nd is a bounded domain with a piecewise smooth boundary r, rid — 2 or 3 represents 
the number of space dimensions, u T — (ui, u 2 , ...u m ) is a vector of m unknown functions of 
x. A,- and A are mxm matrices which depend on x,f is a given vector-valued function, B 
is a boundary operator, and £ is a given vector-valued function on the boundary. Without 
loss of generality we assume that g is a zero vector. 

Here, we do not discuss the existence and uniqueness of the solution to (3.1), because 
these depend on the structure and properties of L and B, and the vector /. In the 
following discussion, it is assumed that the problem (3.1) has a unique solution. We indicate 
that if there is a solution to (3.1), then the following least-squares method produces an 
approximate solution(see,e.g.,[30,31]). 

Throughout, L 2 ( fl) denotes the space of square-integrable functions defined on Cl with 
inner product 

(u,v) = / uvdCl u,v € L 2 (Cl) (3.3) 

in 

and norm 

IMIo = ( u > u ) u€L 2 (n) (3.4) 

We define the Sobolev space as: 

jy^n) = {u e L 2 {ci),d a u € l 2 (ci), v|«| < 1} (3.5) 


where a = (a l5 a 2 > •••> a n d ) € N nd and |a| = a! + a 2 + ••• + <Xn d -> an< ^ define their associated 
norms by 

INI? = E Il a ““ll5 (3-6) 

I«I<1 

For the vector- valued function u with m components, we have the product spaces 

L 2 (n) = (L 2 (n)) m (3.7) 

H}{Cl) = (tf^fi))" 1 (3.8) 

and the corresponding norm 

m 

Hallo = E Kilo (3-9) 

i=i 

m 

Halil = EKIli (3.10) 

j=i 


Considering the boundary condition of the boundary-value problem, we also define 
the function space 

S = {ue{H 1 {Cl)) m ]Bu = 0 onf} (3.11) 
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Let us suppose that / G L 2 and L : 5 — ► L 2 . For an arbitrary trial function u G S, we 
define the residual function R = Lu — /. The LSFEM is based on minimizing the residual 
function in a least-squares sense. 

We construct the least-squares functional 

J(u) = ||Lu - /||o = (L« - /, Lu - /) (3.12) 

Taking variation of I with respect to u, and setting 61 = 0 and 6u = w, lead to the 
least-squares weak statement: Find u G S_ such that 

(Lu>,Lu) = (Lu;,/) Vtu G 5 (3.13) 


In the approximate analysis, we first discretize the domain as a union of finite elements 
and then introduce an appropriate finite element basis. Let N e denote the number of nodes 
for one element and V’y denote the element shape functions. If equal-order interpolations 
are employed, that is, for all unknown variables the same finite element is used, we can 
write the expansion 


Ne 




j = i 




V u m / 


(3.14) 


where (ui,ii 2 , ...,u m )y are the nodal values at the jth node, and h denotes the mesh 
parameter. 


Introducing the finite element approximation defined in (3.14) into the weak statement 
(3.13), we have the linear algebraic equations 


KU = F (3.15) 

where the U is the global vector of nodal values. The global matrix K is assembled from 
the element matrices 

K e = [ (LV>1 , LV>2 > •••> LV’jVe) T (LV’l , L^2 > •••> Ij1pNe)dft (3.16) 

JO' 


in which fl e C fl is the domain of the eth element, and T denotes the transpose, and the 
vector F is assembled from the element vectors 



in which 

L rpj - ipj, x Ai + V’y.yAa + ^>A 3 + V’yA (3.18) 

We observe that the matrix K is symmetric and positive definite. Therefore, the itera- 
tive methods, such as the preconditioned element-by-element conjugate gradient method[27], 
can be employed. This property means that large-scale problems can be effectively solved 
on vector and parallel computers by using LSFEM. 

Here it is important to emphasize that there is no added weighting parameter in our 
LSFEM. Each equation in the first-order system is equally treated. Usually, the govern- 
ing equations are nondimensionlized by using the characteristic quantities of the problem. 
The nondimensionlization or scaling is equivalent to weighting the residuals in the least- 
squares functional. Although the choice of the characteristic quantities is arbitrary, the ex- 
act(physical) solution at which the least-squares functional takes the minimum value(zero) 
is independent of the scaling. In other words, the weighting or scaling can not change the 
minimum point. However, the weighting can change the condition number of the result- 
ing algebraic system obtained from discretization. For the Cauchy-Riemann equations, 
we did both numerical experiments and a theoretical analysis to investigate the effect of 
the weighting factor [25, 26]. It was observed that, when the weighting factor was wildly 
changed from 10“ 8 to 10 8 , the numerical solutions were almost fixed, except the number 
of conjugate gradient iterations was different. When the wieghting factor was equal to one, 
the number of iterations to convergence was minimized. 

The formulation of LSFEM is simple and systematic. Once a LSFEM based on a 
first-order system is coded, adaptation of the program to other problem is simply carried 
out by modifying the subroutines associated with the coefficient matrices Ai, A 2 , A 3 and 
A and the vector function /. In this way, one may develop a general-purpose program. 


4. Incompressible Viscous Flow 

In the past two decades numerous finite element schemes have arisen for solving the 
Navier-Stokes equations of incompressible viscous flow(see,e.g.,[5]). To date, most of these 
scheme require the use of non-equal-order interpolations and the solution of linear system 
of equations with non-symmetric and non-positive-definite matrix. In this section, we 
discuss an application of LSFEM which can circumvent these difficulties. 


4.1 Stokes Problem 

Consider the Stokes problem: Find the velocity u = (tti, 1 * 25 ^ 3 ) and the pressure p 
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such that 

in ft (4.1a) 

in ft (4-16) 

u = 0 on T (4.1c) 

where all variables are nondimensionalized, and Re is the Reynolds number. 

The most popular method for the Stokes problem is based on the Galerkin mixed 
formulation(see,e.g.,[4,5]): Find {u,p} € (.Ho(ft)) 3 x (L 2 (fi)/9t) such that 

(V.u, 9)=0 v g eL 2 (n) 

a(u,v) - (V-v,p) = (f,v) VvG (IfZ(n)) 3 

where 

^(ft) = {u<Eff 1 (n);u = Oon T} 

o(u,u) = v ”») 

*=1 

The weak statement (4.2) can be obtained directly from the stationary condition for 
the Lagrangian 

L{u, p) = ^a(u, u) - (/, u) - ( V • u, p) (4.3) 

The solution {u, p} defines a saddle-point of this functional. 

The existence of finite element approximate solution to (4.2) depends on choosing a 
pair of space V_ h C (Jfo(ft)) 3 and Qh C L 2 (Vl) such that the Ladyzhenskaya-BabuSka- 
Brezzi(LBB) “inf-sup” condition(see [32-34,4]) 

m/ > a > o (4.4) 

qeQh vev, || v ||i|| 9 || 0 

where a is independent of the mesh size h , holds. This condition precludes the use of 
equal-order interpolations. Although for two-dimensional problems quite a few convergent 
pairs of velocity and pressure elements have been developed(see,e.g.,[4,5,34]), most of these 
combinations employ some basis functions that are inconvenient to be implemented. For 
three-dimensional problems, this difficulty becomes more severe and only rather elaborate 
constructions can pass the LBB test. The other basic difficulty associated with the mixed 
method is the non-positiveness of the linear algebraic systems. 


(4.2a) 

(4.26) 

(4.2c) 

(4.2(f) 


V • u = 0 


__ A „ + V p=/ 
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features of the LSFEM are more clearly explained by a simple one-dimensional two-point 
boundary-value problem in [41]. 

In the following, we describle another LSFEM developed by Jiang and Chang[4l]. 
Introducing an auxiliary variable, the vorticity w = V x u, the Stokes equations can be 
written as 

V • u = 0 

-^-Vxw + Vp = / in n (4.7) 

Re ~~ 

oj — V x u = 0 

We shall consider the two-dimensional problem only: 

du dv 
dx dy 


dp 

dy 


dp 1 dw 
dx ' Re dy 


1 

Ye~di = fy 


in 


n 


du dv 

u + ai~ai = 0 


(4.8) 


where ( f x ,fy ) is the body force components, and 0 is a bounded domain in JR 2 with a 
piecewise smooth boundary T. Let (Fi, T 2 , T 3 , I^, Ts) denote the sides of T. The unit 
outward normal vector to T is denoted by n, and the tangential vector to T by t. 


We now write (4.8) in a general form of first-order system (3.1a) , in which 



f 1 

0 

0 

0 \ 



1 

0 

o \ 


0 

0 

1 

0 


0 

0 

0 

1 

Ai = 

° 

0 

0 

1 

Re 

A 2 — 

0 

0 

1 

Re 

0 


VO 

-1 

0 

0 J 


VI 

0 

0 

0 J 


/0 0 0 0 \ 
0 0 0 0 I 
0 0 0 0 I 

VO 0 0 1/ 




(4.9) 


The boundary conditions should be supplemented to complete the definition of the 
boundary-value problem. We may consider the following boundary conditions: 


(a)u, v given on Ti 
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(b) u n ,w given on r 2 

( c) p,ut given on T 3 


(4.10) 


(d) p,tv given on T 4 

( e) p,u n given on 

For example, Ti could be the inflow, outflow or wall boundary; r 2 ,r 3 ,r 4 , and T 5 could 
be the free surface, inflow, or outflow boundary. 

The solvability of this elliptic boundary value problem depends on the combination 
of the boundary conditions. The system of equations is of the Petrovsky type. However, 
for many practical problems, the boundary operator does not satisfy the Lopatinski condi- 
tion[l9]. Thus, the proof of solvability is not trivial. This problem and a theoretical error 
analysis are addressed in [42]. 

Now we are ready to use LSFEM as described in Section 3.1. As usual we use the 
Gaussian quadrature technique to evaluate the coefficients of K e and F e . The least- 
squares method with Gaussian quadrature is equivalent to the collocation least-squares 
method[43]. Therefore, the total number of collocation points(‘Gauss’ points) should be 
compatible with that of unkowns to get good results. For this reason, in our numerical 
experiments we use reduced integration[14]. 

We chose a model problem studied by Oden and Jacquotte[44]. This model problem 
corresponds to the polynomial divergence-free velocity-pressure-vorticity field: 

u(x,y) = x 2 (l - x) 2 (2y - 6 y 2 + 4y 3 ) 

v(x,y) = y 2 ( 1 — y) 2 (— 2x + 6 x 2 -4x 3 ) 

p(x,y)=x 2 -y 2 (4.11) 

w(x,y) = — x 2 (l — x) 2 (2 — 12y + 12y 2 ) + y 2 (l — y ) 2 (— 2 + 12 x — 12 x 2 ) 
with the boundary conditions shown in Figure 4.1. 

We have tested the bilinear element and the eight-node quadratic element using uni- 
form meshes. The one-point Gaussian quadrature for the bilinear element and the 2 x 2 
quadrature for the eight-node quadratic element are employed. The numerical results of 
the rate of convergence are shown in Figure 4.2. It is found that in all tested cases 

(ll u - u h||o + ll v - «fc||o)* < ch k+1 , ||p — Pfcllo < ch k+1 , ||u - Wfc||o < ch k+1 (4.12) 
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where h is the mesh size; as a subscript , h denotes the approximate solution; k is the order 
of the polynomial. That is, all variables u,v,p,u converge in L 2 norm at the optimal rate, 
as predicted by a theoretical error estimate in [42]. 


4.2 Incompressible Navier-Stokes Problem 


The method presented in Section 4.1 can be generalized to solve the two-dimensional 
incompressible Navier-Stokes equations written in the following first-order quasi-linear 
velocity-pressure-vorticity formulation[45]: 

du dv 
dx dy 


du du dp 1 du 
dx V dy + dx + Re dy x 
dv dv dp 1 du _ 
dx+ V dy dy Re dx v 
du dv 

u + — — — = 0 

dy dx 


in 0 


(4.13) 


The boundary conditions should be supplemented to complete the boundary-value prob- 
lem. 

We now write (4.13) in a general form of first-order system (3.1a) , in which 



1 

0 

0 

0 



(0 

1 

0 

0 

\ 

Ai = 

u 

0 

0 

u 

1 

0 

0 

1 

Re 


A 2 = 

V 

0 

0 

V 

0 

1 

1 

Re 

0 



^0 

-1 

0 

0 

) 


U 

0 

0 

0 

) 


A = 


0 0 0 0 ' 

0 0 0 0 

0 0 0 0 

.0 0 0 1 


\ 



/u\ 

f = 

/* ’ 

U = 

V 


f« j 


p 

f 

0 ) 


\w/ 


(4.14) 


A successive substitution method is used to solve this quasi-linear problem. The 
corresponding solution of the Stokes problem is taken as the initial guess. The velocity 
field at the previous step is used to calculate the coefficient matrices. 

The method is tested for the driven cavity problem by using 40 x 40 nonuniform 
bilinear elements. The difference between the results of two successive steps is defined as 


e = 


max 

i—l,N n 


u" — < +1 
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where * denotes the node, N n is the total number of nodes, and n denodes the substitution 
level. The substitution continues until the difference e becomes less than the tolerance 
10 -4 . The required number of iterations are 8, 12, 15, 57, and 160 for the Reynolds 
number Re of 100, 400, 1000, 3200, and 5000, respectively. The numerical results compare 
favorably with a finite difference solution using fine grid. Figure 4.3 illustrates the mesh 
used and Figure 4.4 compares the LSFEM results with those obtained by Ghia[46], Figure 
4.5 through 4.8 show the flow pattern, the streamlines, the pressure contours and the 
vorticity contours obtained with the LSFEM for Re = 5000. 

We conclude this section by comparing LSFEM with other methods. In contrast to the 
Galerkin mixed method based on the velocity- pressure formulation[5,34], LSFEM based 
on velocity-pressure- vorticity formulation does not depend on the LBB condition , can 
accommodate equal-order interpolations, and always produces a symmetric, positive and 
definite matrix. In contrast to the penalty method[5], LSFEM does not have any added 
parameter in the scheme. In other words, LSFEM is robust. In contrast to the stream 
function- vorticity formulation [5, 34], LSFEM produces the velocity and pressure directly 
and can be extended to three-dimensional problems. 


5. Transport 

In this section, we consider the two- or three-dimensional convection (hyperbolic) prob- 
lems: 


d<t> 

-y + u • V4> 
at 

= 0 in n x (0, T ) 

(5.1a) 

4> = g 

on T_ x (0, T) 

(5.16) 

o 

II 

in Cl for t = 0 

(5.1c) 


where <f> is, for example, the concentration, the convective velocity vector u is assumed 
known and divergence-free in H, and the boundary data g is prescribed only on the inflow 
boundary T_ = {x € r;n(x) • u(x) < 0} in which n(x) is the outward unit normal to T at 

x G r. 

The classic Galerkin variational problem is to find <j> G S — (JT 1 (n); t > — g on T_} for 
each t, such that 

K ^ +u- V<£) = 0 Vt l )6V = {fl 1 (fl);t; = 0onr.} (5.2) 

(j I 

together with the initial condition 

HO) = <f> o 
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For convective transport problems!, the classic Galerkin finite element method suffers 
from spurious oscillations in the same way as the central finite differencing. Motivated 
by backward differencing to reflect the ‘upstream’ dependence of the solution, the Petrov- 
Galerkin method and the streamline upwind Petrov-Galerkin(SUPG) method were devel- 
oped(see,e.g.,[47,8]). The weak statement of the SUPG method is to find <)> € S for each 
t, such that 

(w + tu-Vw, — + u-V<^)= 0 Vto € V (5.3) 

where r is a parameter. The term “Petrov-Galerkin” is used to indicate that the test 
functions are no longer the same as the trial functions for the approximation. The term 
“streamline upwind” comes form the term u • Vto in the test functions. 

In the steady case, the SUPG formulation is: Find <f> 6 S such that 

(tu + ru • Vtu, u • V<j>) = 0 Vtu (E V (5.4) 

Note that this formulation combines the usual Galerkin term and the least-squares term 
together. For this reason, Hughes has renamed the SUPG method as the Galerkin/least- 
squares method[48]. By fine-tuning the upstream weighting parameter r in the scheme, 
excellent results can be obtained. 

Recently, the ideas of the Lax-Wendroff scheme in finite difference method have mo- 
tivated development of the related Taylor-Galerkin finite element method[49,6]. Let us 
briefly describe the Taylor-Galerkin method. Since equation (5.1a) defines an evolution 
statement, there must exist the Taylor expansion 

— <j> n + A t<j>™ + —A t* <t>tt “b ••• (5*5) 

where the subscipt H ’ denotes the order of temporal derivative at t n , and < n +i = t n + At. 
Equation (5.1a) permits restatement of the first derivative term and the second derivative 
term in (5.5) as 

fa = -« • V<£ (5.6) 

fa t = -u-Vfa = u-V{u-V<f>) (5.7) 

Substituting (5.6) and (5.7) into (5.5) yields 

<f> n+1 =<f> n - A tu ■ Vfa 1 + \ At 2 u • V(u • V<T) + (5.8) 

Testing (5.8) against w € V and integrating by parts we obtain the Taylor-Galerkin weak 
statement: Find <f> € S such that 

(w, <t> n+l ) = (w, <f> n — A tu • V<j> n ) — ^ At 2 (u • Vtu, u ■ V^ n ) 
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Vu; G V 


(5.9) 


+ ^A t 2 (w, (u • V<£ n )(u • n)) r+ + ... 

where T + denotes the outflow boundary. 

Since the term (tu, <j> n+1 ) introduces the mass matrix, the Taylor-Galerkin method 

(5.9) requires inversion of the mass matrix. To maintain the explicit character of the 
method, one may use the lumped (diagonal) mass matrix. The method thus suffers from 
more severe stability limitations(Courant number C =| u \ At/h < l/\/3). An implicit 
Taylor-Galerkin method[50] can be constructed, but it requires solution of nonsymmetric 
matrix. Furthermore, the Taylor-Galerkin method is essentially a high-order scheme, and 
hence may promote oscillations at discontinuities. Artificial dissipation can be introduced 
to suppress oscillations, but the form of this added term is not unique and the associated 
parameters must be specified. 

For first-order hyperbolic systems, another choice is LSFEM[51-54]. This method uses 
the backward scheme or the Crank-Nicolson scheme of finite differences in the temporal 
domain and finite elements in the spatial domain. For one-dimensional problem, it was 
demonstrated that LSFEM is unconitionally stable for all Courant numbers, and naturallly 
acts in a manner similar to upwinding. 

Here let us apply LSFEM to the problem (5.1). Backward differencing (5.1) through 
a time step and writing <j> n for <j>(x,t n ), the resulting semidiscrete problem is 

<j> n+1 _ + At ^ . v<£ n+1 = 0 (5.10) 

Now we are ready to use LSFEM as described in Section 3. In order to explain how the 
method works, we deduce the least-squares weak formulation. 

For any admissible 4> n+x {x ) 6 S , we obtain from (5.10) the residual R, and the least- 
squares functional I = J n R 2 dft becomes 

I(^ n+1 ) = ||<T +1 - <t> n + A tu • V<f> n+X 1|2 (5.11) 

Applying the stationary condition SI = 0, and setting 6<j> n+1 = w, we obtain the weak 
statement for (5.10): Find <f> n+1 £ S such that 

( w + A tu • Vw, <f> n+1 + A tu • V<£ n+1 ) = (u; + A tu • Vtu, <f> n ) Vu; € V (5.12) 

Clearly (5.12) can be interpreted as a Petrov-Galerkin weighted-residual statement for 

(5.10) with test function w — (1 + Atu- V)w. More importantly, we note that the convective 
velocity u enters explicitly in this test function so that this form is analogous to SUPG 
methods. Hence, the effect of the least-squares procedure in two or three- dimentions is to 
add numerical dissipation preferentially in the local flow direction in a manner similar to 
SUPG method. 
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Introducing the increment A <f> = 4> n+l — <f> n , the formulation (5.12) may be rewritten 
in the A-form: Find A <f> £ V such that 

(w + Atu-Vw, A<f> + Atu • V(A0)) = — Ai(iu + Atu • Vu>, u • V0 n ) Ww£V (5.13) 

The formulation (5.13) clearly shows that the resulting matrix of LSFEM is symmetric. 

If the boundary condition (5.1b) is independent of the time t. The steady state solution 
exists. In this case, when the steady state is reached by the time marching, A <f> becomes 
zero, and the corresponding weak statement is: Find <j> 6 S such that 


[w + Atu • Vw, u • V0) = 0 Vtu 6 V 


(5.14) 


Note that the formulation (5.14) is the same as the SUPG or Galerkin/least-squares for- 
mulation (5.4). Although the Galerkin/least-squares formulation can be deduced by the 
least-squares procedure, there is a significant difference between these two methods. Even 
for steady state cases, LSFEM uses (5.12) or (5.13) as a time marching procedure. For each 
time-step, the algebraic system obtained by LSFEM is symmetric, while the Galerkin/least- 
squares method leads to a nonsymmetric matrix. 


We should also point out that the success of LSFEM for the hyperbolic systems 
depends on the choice of time step At. If the time step is too large, the numerical results at 
discontinuities will be diffused too much; if the time step is too small, the solution will have 
significant oscillation at discontinuities. Our numerical experiments tells us that one may 
choose the time step such that the Courant number C ~ 10 ~ 50. Of course, a theoretical 
investigation about this problem is needed. Although LSFEM has this disadvantage, the 
symmetry of the matrix makes the method attractive. 


As a numerical example, we consider the water-oil displacement problem in 5-spot 
pattern. A displacement front propagates from the water injection well at the lower- 
left corner towards the oil production well at the top-right corner. The flow domain is 
assumed homogeneous, and the velocity field is then given by the solution of Laplace’s 
equation for this source-sink system with symmetry (no flow) conditions on the sides to 
enforce periodicity as 

x x — 1 

x 2 + y 2 ( x - l) 2 + (y - l) 2 


y y - i 

x 2 +y 2 (x - l) 2 + (y - l) 2 


(5.15) 


Using this specified velocity field in the least-squares finite element analysis, the finite- 
element system follows directly from the introduction of the finite element expansion for 
<f> into (5.12). The unit square domain is discretized to a uniform 40x40 grid of bilinear 
elements with the initial data 0(x,y, 0) = 0 and the injection value 0(0,0, t) = 1 for t > 0. 
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The equispaced solution contours for concentration <j>(x,y,t) are shown in Figure 5.1 at 
t = 0.05 and 0.3 and the fixed time step At = 0.0025. We note that there is numerical 
dissipation of the propagating front and no oscillation. 


6. High-Speed Compressible Flow 

Recently, much attention has been focused on the development of finite element al- 
gorithms for the analysis of high speed compressible flows especially for treating discon- 
tinuities. Although considerable success has been achieved by using the methods such as 
the streamline-upwind Petrov-Galerkin method[8], the Taylor-Galerkin method[55,7], the 
Taylor-Galerkin method with flux-corrected transpose(FCT) and mesh refinement [56-59], 
the block relaxation via Godunov’s method[60] ,the characteristic Galerkin method[6l], and 
the non-osillatory discontinuous Galerkin method [62], more investigation is still needed to 
compare the efficiency and accuracy of these methods and alternative schemes. 

For the numerical solution of the one- and two-dimensional compressible Euler equa- 
tions, we proposed a class of least-squares methods[63,64|. We begin by considering 
the first-order implicit time-differenced non-conservative formulation. The least-squares 
method is then employed to minimize the residual in the L 2 or H 1 norm. It was demon- 
strated that the associated numerical dissipation arises naturally as a consequence of the 
method. 


6.1 L 2 Method 

Now let us describe LSFEM for the two-dimensional compressible Euler equations in 
the form of a first-order system: 


du du 

~dt + Al dx 



0 


where u = (p, u,v,p), and 


(u 

p 

0 

0 
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( v 

0 

p 
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II 
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0 
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IP 

0 
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\o 

0 

IP 

V 
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( 6 . 1 ) 


in which p is the density, (u,v) are the fluid velocity components , p is the pressure, and 
7 is the specific heat ratio. The boundary conitions and the initial condition are needed 
to complete the description of the problem. 
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For a given time step At = f n+1 — t n , we linearize the problem by setting A” = 
Ai(u n ),A 2 = A 2 (u"). Backward differencing leads to the implicit time-differenced prob- 
lem 


u n+1 _ u n + AtA5*^| + AtA$ ?= + . = o 

dx 


dy 


( 6 . 2 ) 


Now we are ready to use the formulation of LSFEM in Section 3 to get the solu- 
tion u n+1 for each time-step. Here, the first-order partial differential operator L has the 
following form: 

L = At A" + At A£ -J- + E (6.3) 

ox ay 

where E is an identity matrix. 

We proceed further to discern why this method can capture shocks. For this reason 
we deduce the weak statement. 

The basic least-squares method for the system (6.2) amounts to minimizing the L 2 
norm of the residual R for admissible u n+1 in (6.2); i.e., minimizing the objective functional 

Io = WRWl (6.4) 


where 


R = u n+1 -u n + AfA?%^ + A<A£ — 

ox 


dy 


Taking variations of Jo with respect to u ,l+1 and setting uj = 6u n+1 , 6Io = 0 leads to the 
least-squares weak statement 

<* + «*£ + A ‘ A # ^ + A ‘ A ?^ + A ‘ A ^> 


= {tv + AtA" + AfA£^, u n ) 
ox dy 


(6.5) 


When the steady state is reached by time marching, the corresponding weak statement 
becomes 


/ . dw . dw. . . . dw dw du . du , „ .. 

+ A2 ^ ) + A ‘ (Al fe +A2 ai7’ Ai ^ +A2 ^ )_q (6 ' 6) 


Remark For steady state problems we still use the formulation (6.5) as the time 
marching procedure, since the resulting algebraic system is symmetric. The formulation 
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(6.6) is not appropriate for practical calculation, because it is a nonlinear problem and the 
resulting matrix is nonsymmetric. 

Let us examine the structure of (6.6). The first term in (6.6) corresponds to the 
standard Galerkin weak statement. The second term leads to a symmetric, positive and 
definite matrix, and always acts as a numerical dissipation to stabilize the soluton by 
smoothing out any discontinuities. 

We now demonstrate some features of the L 2 LSFEM using three numerical examples. 

A standard test problem corresponding to the reflection of a shock from a wall is 
depicted in Figure 6.1. On the upper boundary of the flow domain p = 1.7, u = 2.6185, v = 
— 0.5082, p = 1.5282 and on the upstream boundary p — 1.0, u = 2.9, v = 0.0, p = 0.7143 
are prescribed. The lower boundary, where v = 0, is the wall from which the generated 
shock reflects, and the downstream boundary conditions remain free for outflow. The 
initial data were prescribed as constant at values given on the upper boundary and the 
specific heat ratio is 7 = 1.40. In the calculation, a uniform 20x60 mesh of bilinear 
elements was used. The solution was integrated with time-step At = 0.33333 until an 
essentially steady state is obtained in 12 time-steps. The pressure contours for this steady 
solution are given in Figure 6.1. Qualitatively, it is seen that the flow physics are correctly 
modelled. Although the shock is somewhat smeared, the oscillations are absent and the 
calculation is efficient. 

The second problem is a Mach 3 flow (with 7 = 1.40) over a 20° ramp. The gas enters 
with uniform flow conditions through the left boundary of the domain and an oblique shock 
develops at the root of the ramp. A mesh with 824 bilinear elements and the computed 
pressure contours are illustrated in Figure 6.2. In the calculation, the initial data were 
prescribed as constant at the value given on the left boundary, and the time-step was 
At = 0.33333. The steady state was obtained in 16 time-steps. 

We also considered a cylinder in a supersonic flow with Mach number M 0 0 = 2, 7 = 
1.40. Figure 6.3 presents the mesh with 800 bilinear elements and the computed pressure 
contours. 

We note that in all three numerical examples , LSFEM with bilinear elements produces 
non-oscillatory shock profiles and the results compare favorably with those in the literature 
based on other finite element methods. 


6.2 H l Method 

Numerical experiments show that as long as the time-step At is large enough ( the 
corresponding Courant number ~ 10 - 50 ), the L 2 method with linear elements gives 
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non-oscillatory shock profiles in the computed steady state. However, for high-order ele- 
ments, computed solutions based on this method may have oscillations. The behaviour of 
this method is similar to that of high-order finite difference schemes. To circumvent this 
problem, one may use low-order (linear) elements around discontinuities and high-order 
elements elsewhere. However, the implementation is not so easy. Thus, we proposed an 
alternative technique in which the objective function is modified to control the residual 
derivatives also[63-65|. Accordingly, the following multi-objective optimization problem is 
constructed [66]: 



(6.7) 


in which the minimization of Jo is the main objective. As a technique to solve (6.7) , we 
minimize the modified functional 


/=ll£HS + /?/* 2 [ll^ll5 + ll^ll2] (6.8) 

where /? is a small parameter, 0 < (3 C 1. Note that minimizing the functional (6.8) can be 
interpreted as minimizing the weighted H l norm of R. However, the form in (6.8) contains 
second derivatives of u, which implies that, for conformity, C 1 elements are appropriate. 
Generally, for two- and three-dimensional problems, simpler elements are desirable. An 
approximate formulation based on (6.8) can be constructed using high order C° elements, 
we may approximate (6.8) as 


T t T~» 7-.\ /dR tdR\ 1 

h-{R, R)+/3h [(— , + (^)nl ( 6 * 9 ) 

where Cl denotes element interiors. For example, with C° quadratic or cubic elements, the 
contributions of the second derivatives in element interiors can be calculated for (6.9). 

We now demonstrate the shock capturing ability of the H 1 LSFEM using the same 
shock wave reflection problem as in Section 6.1. The computed pressure contours with 
8-node quadratic elements are shown in Figure 6.4. Comparing Figure 6.4 with Figure 6.1, 
we observe that for the same grids this new techniqure gives the better resolution of the 
shocks. 


6.3 Conservative Form 

The LSFEM described in Section 6.1 is based on the nonconservative formulation. 
Here we would like to construct a conservative LSFEM (in the sense of the steady state). 
The Euler equations governing two-dimensional compressible inviscid flows can be written 
in conservative form as 
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( 6 . 11 ) 


where 


dg dE_ dQ 
dt dx dy ~ 



p \ 


/ pxi \ 

pu 

F = 

pu 2 + p 

pv 


puv 

pe J 

(pe -f p)uj 


G = 


pv 

puv 

pv 2 +p 
V ( pe + p)v 


( 6 . 12 ) 


in which e is the total energy, and for the case of perfect gas the equation of state is 


P= ft - 1)(= - + V 2 )) 


(6.13) 


Retaining the conservative variables as an intermediate step, we first convert ( 6 . 11 ) 
(temporarily) into the following nonconservative form: 


dq — dq 
dt dx 

- dq 

+ A 2 — — 0 
dy 



(6.14) 

( 0 

1 

0 



j \l(u 2 + v 2 ) - u 2 
I —uv 

(3 -7)« 
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:) 

(6.15) 

V {l{u 2 + v 2 ) — 7e)u 6 — 7 u 2 

— 7Ut; 

7 u J 



A2 


° 

—uv 

£71 (u 2 + v 2 ) - v 2 
V (7(u 2 + v 2 ) - 'ye)v 


0 1 0 \ 

v u 0 I 

-7*» (3 - l)v 7 I 

— ')uv l — 7V 2 7 v J 


in which 


7 = 7-1, 


c = 7e - 7 


(u 2 + v 2 ) 


(6.16) 


Using the same procedure as before we obtain the coresponding least-squares weak 
statement similar to (6.5). Introducing A q = q n+1 — q n and regrouping the terms, we have 
then a conservative least-squares weak statement: 


, — n dw — n du) — n 3Ao — n <9A<7 

(w + AtA x + AiA 2 + AtAj — — h AtA 2 


A . . — n dw . — n dt£ 

-At(w + AtA x (- At A 2 — , 

ox oy 


dF n 


dx 


dy 


(6.17) 
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Here A q is the unknown increment in conservative variables for the time step At. We 

may use the conservative variables at the previous time step to calculate the nodal values 

of components of flux F n and G n , then use a finite element approximation to calculate 
dF n d G n 

-gj- and -gp- . Once the steady conservation law 

dE_ dQ _ 
dx dy “ 

is satisfied, the increment A q becomes zero, and the calculation is terminated. Our nu- 
merical experiments on the shock wave reflection problem show that there is no essential 
difference between the results of nonconservative and conservative LSFEM. More numerical 
experiments on different problems are needed to compare these two methods. 



7. Conclusions 

A variety of complex flow problems have been solved using the least-squares finite ele- 
ment method. Perhaps no other single method can solve the same range of flow problems. 
This paper has demonstrated that the method is simple, general, and reliable. It is hoped 
that this review will help stimulate the development of a general-purpose flow solver based 
on the least-squares finite element method , which in turn will speed the development of 
new applications. 
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Figure 4.1 Boundary conditions for Stokes problem 



Figure 4.2 Computed convergence rate for Stokes problem. 



Figure 4.3 Finite element mesh (40x40 bilinear elements) for the cavity flow. 
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Figure 4.4 Horizontal velocity profile for cavity flow at Reynolds number 5000 
present (40x40 grids), A Ghia[47] (257x257 grids). 
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Figure 4.5 Flow pattern of cavity flow at Re=5000 
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Figure 5.1 


(b) at t = 0.30 

Concentration contours for five-spot problem 
(40 x 40 bilinear elements, At = 0.0025) 
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Figure 6.1 Pressure contours for wave reflection problem ( L 2 method, 20x60 bilinear 
elements, At = 0.33333, number of steps = 12 ) 
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Figure 6.3 Finite element mesh and pressure contours for the Mach 2 flow past a circular 
cylinder ( L 2 method, 800 bilinear elements, A t — 0.1, number of steps=36) 



Figure 6.4 Pressure contours for wave reflection problem ( H l method, P = 0.0001, 
10x30 8-node quadratic elements, A t = 0.33333, number of steps = 14 ) 
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